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Abstract 

A  computational  methodology  is  developed  to  efficiently  perform  uncertainty  quan¬ 
tification  for  fluid  transport  in  porous  media  in  the  presence  of  both  stochastic 
permeability  and  multiple  scales.  In  order  to  capture  the  small  scale  heterogeneity, 
a  new  mixed  multiscale  finite  element  method  is  developed  within  the  framework 
of  the  heterogeneous  multiscale  method  (HMM)  in  the  spatial  domain.  This  new 
method  ensures  both  local  and  global  mass  conservation.  Starting  from  a  specified 
covariance  function,  the  stochastic  log-permeability  is  discretized  in  the  stochastic 
space  using  a  truncated  Karhunen-Loeve  expansion  with  several  random  variables. 
Due  to  the  small  correlation  length  of  the  covariance  function,  this  often  results  in  a 
high  stochastic  dimensionality.  Therefore,  a  newly  developed  adaptive  high  dimen¬ 
sional  stochastic  model  representation  technique  (HDMR)  is  used  in  the  stochastic 
space.  This  results  in  a  set  of  low  stochastic  dimensional  subproblems  which  are 
efficiently  solved  using  the  adaptive  sparse  grid  collocation  method  (ASGC).  Nu¬ 
merical  examples  are  presented  for  both  deterministic  and  stochastic  permeability 
to  show  the  accuracy  and  efficiency  of  the  developed  stochastic  multiscale  method. 
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1  Introduction 


Flow  through  porous  media  is  ubiquitous,  occurring  from  large  geological 
scales  down  to  microscopic  scales.  Several  critical  engineering  phenomena  like 
contaminant  spread,  nuclear  waste  disposal  and  oil  recovery  rely  on  accu¬ 
rate  analysis  and  prediction  of  these  multiscale  phenomena.  Such  analysis  is 
complicated  by  heterogeneities  at  various  length  scales  as  well  as  inherent 
uncertainties.  For  these  reasons  in  order  to  predict  the  flow  and  transport 
in  stochastic  porous  media,  some  type  of  stochastic  upscaling  or  coarsening 
is  needed  for  computational  efficiency  by  solving  these  problems  on  a  coarse 
grid.  However,  most  of  the  existing  multiscale  methods  are  realization  based, 
i.e.  they  can  only  solve  a  deterministic  problem  for  a  single  realization  of  the 
stochastic  permeability  field.  This  is  not  snfficient  for  uncertainty  quantifica¬ 
tion  since  we  are  mostly  interested  in  the  statistics  of  the  flow  behavior,  snch  as 
mean  and  standard  deviation.  In  this  paper,  we  propose  a  stochastic  multiscale 
approach  which  resolves  both  uncertainties  and  subgrid  scales  by  developing  a 
new  multiscale  method  and  adopting  a  newly  developed  adaptive  high  dimen¬ 
sional  stochastic  model  representation  technique  (HDMR).  The  goal  of  the 
multiscale  method  is  to  coarsen  the  flow  equations  spatially  whereas  HDMR 
is  used  to  address  the  curse  of  dimensionality  in  high  dimensional  stochastic 
spaces. 

One  of  the  challenging  mathematical  issnes  in  the  analysis  of  transport  through 
heterogeneous  random  media  is  the  multiscale  nature  of  the  property  varia¬ 
tions.  Complete  response  evaluation  involving  full-scale  spatial  and  temporal 
resolution  simulations  of  multiscale  systems  is  extremely  expensive.  Compu¬ 
tational  techniques  have  been  developed  that  solve  for  an  appropriate  coarse- 
scale  problem  that  captnres  the  effect  of  the  subgrid-scales.  The  most  popu¬ 
lar  techniques  developed  for  such  upscaling  fall  under  the  category  of  multi¬ 
scale  methods  viz.  the  multiscale  finite  element  (MsFEM)  method  [1,2],  the 
variational  multiscale  (VMS)  method  [3,4]  and  the  heterogeneous  multiscale 
(HMM)  method  [5,6].  The  MsFEM  was  originally  developed  in  [1,2]  for  the 
solution  of  elliptic  equation  based  problems  with  multiscale  coefficients  using 
conforming  linear  hnite  elements.  The  primal  unknown  is  the  nodal  value, 
e.g.  the  pressure,  and  one  can  obtain  the  velocity  by  calculating  the  gradient 
of  the  pressure  field  given  the  finite  element  solution.  The  result  is  generally 
not  accurate  and  conservation  of  the  flux  in  each  element  may  be  violated, 
which  is  an  important  property  for  the  numerical  solution  of  transport  equa¬ 
tions  in  porous  media.  Therefore,  a  mixed  multiscale  finite  element  method 
(MMsFEM)  that  guarantees  the  local  mass  conservation  at  the  element  level 
was  proposed  in  [7]  using  the  lowest-order  Raviart- Thomas  mixed  finite  ele¬ 
ment  [8] .  The  basic  idea  of  the  method  is  to  construct  the  multiscale  finite  ele¬ 
ment  basis  functions  that  incorporate  the  small  scale  information  through  the 
solution  of  a  local  problem  in  each  element  and  couple  them  through  a  global 
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formulation  of  the  problem.  However,  this  work  only  produces  a  globally  mass 
conserving  velocity  field.  This  work  was  extended  in  a  number  of  important 
ways  to  guarantee  mass  conservation  on  both  fine-  and  coarse-scales  [9,10]. 
A  similar  framework  utilizing  the  finite  volume  method  as  the  global  solver 
was  also  proposed  in  [11-13],  which  also  preserves  mass  conservation  at  both 
scales.  The  basic  idea  of  the  VMS  method  is  to  invoke  a  multiscale  split  of 
the  solution  into  a  coarse- scale  part  and  a  subgrid  component.  The  variational 
coarse-scale  problem  is  performed  and  solved  using  the  solution  of  the  local¬ 
ized  subgrid  problem.  Parallel  to  MMsFEM,  a  mixed  finite  element  version  of 
VMS  was  also  proposed  in  [14-16],  which  is  often  called  “Numerical  subgrid 
upscaling”.  A  thorough  comparison  of  the  above  three  methods  for  elliptic 
problems  in  porous  media  flows  can  be  found  in  [17]. 

Unlike  the  MsFEM  which  was  built  on  the  finite  element  method  (EEM), 
the  HMM  is  a  more  general  methodology  for  multiscale  PDEs  (see  [6]  for  a 
review).  The  basis  idea  of  HMM  consists  of  two  components:  selection  of  a 
macroscopic  solver  and  estimating  the  needed  macroscale  data  by  solving  lo¬ 
cally  the  fine-scale  problem.  It  allows  two  different  sets  of  governing  equations 
on  macro-  and  micro-scales,  e.g.  atomistic  simulation  on  micro-scale  and  con¬ 
tinuum  simulation  on  macro-scale  [18,19].  This  framework  was  utilized  to  solve 
multiscale  elliptic  problems  with  the  conforming  linear  FEM  (FeHMM)  [20- 
22].  The  method  was  analyzed  in  a  series  of  papers  [23-25].  However,  unlike 
the  MMsFEM,  there  is  no  discussion  of  the  mixed  version  of  FeHMM  except 
the  work  in  [26],  where  the  author  first  developed  the  theory  of  the  mixed  fi¬ 
nite  element  version  of  HMM  for  the  elliptic  problem  and  proved  the  stability 
and  convergence  of  this  new  method.  However,  the  primitive  idea  in  [26]  is 
only  a  simple  extension  to  the  original  theory  of  HMM  which  in  general  is  not 
suitable  for  realistic  problems  such  as  flow  through  porous  media.  In  addition, 
no  numerical  implementation  was  given  in  [26] .  Motivated  by  the  work  in  [26] , 
in  this  paper,  we  first  develop  and  implement  the  mixed  finite  element  version 
of  HMM  with  application  to  flow  transport  in  heterogeneous  porous  media, 
which  we  will  call  it  mixed  heterogeneous  multiscale  method  (MxHMM). 

All  of  the  above  mentioned  multiscale  analyses  of  such  systems  inherently 
assume  that  the  complete  multiscale  variation  of  the  permeability  is  known. 
This  assumption  limits  the  applicability  of  these  frameworks  since  it  is  usually 
not  possible  to  experimentally  determine  the  complete  structure  of  the  media 
at  small  scales.  One  way  to  cope  with  this  difficulty  is  to  view  the  permeabil¬ 
ity  variation  as  a  random  field  that  satisfies  certain  statistical  correlations. 
This  naturally  results  in  describing  the  physical  phenomena  using  stochastic 
partial  differential  equations  (SPDEs).  The  development  of  efficient  stochastic 
methods  that  are  applicable  for  flow  in  porous  media  has  drawn  significant 
interest  in  the  last  few  years.  Several  techniques  like  generalized  polynomial 
chaos  expansions  (gPC)  [27-29],  perturbation/moment  equation  methods  [30- 
33]  and  stochastic  collocation  method  [31,34-37]  have  been  considered.  Among 
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these  methods,  the  collocation  methods  share  the  fast  convergence  of  the  gPC 
method  while  having  the  decoupled  nature  of  Monte  Carlo  (MC)  sampling. 
This  framework  represents  the  stochastic  solution  as  a  polynomial  approxima¬ 
tion.  This  interpolant  is  constructed  via  independent  function  calls  to  the  de¬ 
terministic  problem  solver  at  different  interpolation  points  which  are  selected 
based  on  special  rules.  Choice  of  collocation  points  include  tensor  product  of 
zeros  of  orthogonal  polynomials  [34,38]  or  sparse  grid  approximations  [39-41]. 
It  is  well  known  that  the  global  polynomial  interpolation  cannot  resolve  lo¬ 
cal  discontinuity  in  the  stochastic  space.  Its  convergence  rate  still  exhibits  a 
logarithmic  dependence  on  the  dimension.  For  high-dimensional  problems,  a 
higher-interpolation  level  is  required  to  achieve  a  satisfactory  accuracy.  How¬ 
ever,  at  the  same  time,  the  number  of  collocation  points  required  increases 
exponentially  for  high-dimensional  problems  (>  10).  Therefore,  its  computa¬ 
tional  cost  becomes  quickly  intractable.  This  method  is  still  limited  to  a  mod¬ 
erate  number  of  random  variables  (5  —  10).  To  this  end,  Ma  and  Zabaras  [42] 
extended  this  methodology  to  adaptive  sparse  grid  collocation  (ASGC).  This 
method  utilizes  local  linear  interpolation  and  uses  the  magnitude  of  the  hier¬ 
archical  surplus  as  an  error  indicator  to  detect  the  non-smooth  region  in  the 
stochastic  space  and  thus  place  automatically  more  points  around  this  region. 
This  approach  results  in  further  computational  gains  and  guarantees  that  a 
user-defined  error  threshold  is  met.  However,  this  method  is  still  not  suit¬ 
able  for  heterogeneous  porous  media  with  small  correlation  length  leading  to 
high  stochastic  dimensionality.  In  recent  work,  Ma  and  Zabaras  [43]  combined 
the  ASGC  with  the  adaptive  stochastic  high  dimensional  model  representa¬ 
tion  (HDMR)  technique  [44] .  HDMR  represents  the  model  outputs  as  a  finite 
hierarchical  correlated  function  expansion  in  terms  of  the  stochastic  inputs 
starting  from  lower-order  to  higher-order  component  functions.  HDMR  is  ef¬ 
ficient  at  capturing  the  high-dimensional  input-output  relationship  such  that 
the  behavior  for  many  physical  systems  can  be  modeled  to  a  good  accuracy 
only  by  the  first  few  lower-order  terms.  An  adaptive  version  of  HDMR  is  also 
developed  to  automatically  detect  the  important  dimensions  and  construct 
higher-order  terms  using  only  the  important  dimensions.  The  heterogeneity  of 
the  porous  media  is  often  due  to  the  small  correlation  length  of  the  covariance 
structure.  All  the  above  mentioned  works  did  not  take  into  account  the  mul¬ 
tiscale  nature  of  the  permeability  Therefore,  in  this  work,  we  will  use  both  of 
these  developments  in  the  stochastic  space  together  with  the  newly  developed 
MxHMM  for  the  spatial  discretization. 

There  exists  several  new  stochastic  multiscale  methods  for  elliptic  problems. 
In  [45]  and  [46] ,  the  variational  multiscale  method  was  extended  to  a  stochas¬ 
tic  version  using  gPC  and  stochastic  collocation  method  respectively  to  solve 
a  simple  diffusion  problem.  The  stochastic  multiscale  finite  element  was  also 
developed  in  [47]  however  only  an  elliptic  problem  was  solved  to  find  the  hy¬ 
draulic  head.  More  related  work  can  be  found  in  [48-50].  In  [48],  the  stochastic 
numerical  subgrid  upscaling  method  was  also  developed  for  the  solution  of  the 
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mixed  form  of  the  Darcy’s  equation  using  the  stochastic  collocation  method. 
However,  in  that  work,  only  the  statistics  of  the  coarse-scale  velocity  and  pres¬ 
sure  were  solved  and  no  flow  transport  problem  was  investigated.  In  [49],  a 
projection  method  for  the  solution  of  stochastic  mixed  multiscale  finite  ele¬ 
ment  method  was  introduced  where  the  velocity  solution  was  projected  onto 
a  multiscale  velocity  basis  functions  which  are  precomputed  from  an  arbi¬ 
trary  number  of  random  realizations.  It  generally  involves  the  solution  of  a 
large  linear  system  of  equations  to  find  the  projection  coefficients  if  the  num¬ 
ber  of  realizations  is  large.  For  each  new  permeability  sample,  this  method 
needs  to  solve  one  coarse-scale  problem  again  and  is  generally  computation¬ 
ally  expensive.  In  addition,  this  method  cannot  provide  the  statistics  of  the 
saturation  directly  and  thus  this  information  was  not  reported  in  their  work. 
In  [50],  this  framework  was  used  to  sample  the  permeability  given  measure¬ 
ments  within  the  Markov  chain  Monte  Carlo  method  (MCMC)  framework  and 
again  no  statistics  of  the  saturation  were  reported.  However,  in  real  reservoir 
engineering,  we  are  primarily  interested  in  mean  behavior  and  a  measure  of 
uncertainty,  e.g.  standard  deviation,  in  the  saturation  of  each  phase.  By  using 
the  adaptive  HDMR  and  ASGC  developed  in  [43],  we  can  obtain  not  only  a 
surrogate  model  for  the  saturation  profile  but  also  can  extract  the  statistics 
of  the  saturation  easily.  Therefore,  the  novelle  contributions  of  this  paper  are 
as  follows:  (1)  We  develop  a  new  mixed  finite  element  version  of  the  heteroge¬ 
neous  multiscale  method  for  the  simulation  of  flow  through  porous  media  in 
the  spatial  domain;  (2)  We  utilize  the  newly  developed  HDMR  technique  to 
address  the  curse  of  dimensionality  that  occurs  naturally  in  this  problem  due 
to  the  heterogeneity  of  the  permeability;  (3)  Finally,  we  investigate  the  effect 
of  the  stochastic  permeability  on  various  statistics  of  the  saturation  using  the 
recently  developed  adaptive  HDMR  method. 

This  paper  is  organized  as  follows:  In  the  next  section,  the  mathematical  frame¬ 
work  of  stochastic  porous  media  flow  problem  in  the  mixed  form  is  considered. 
In  Section  3.2,  the  ASGC  and  HDMR  methods  for  solving  SPDEs  are  briefly 
reviewed.  In  Section  4,  the  theory  of  MxHMM  is  developed.  Various  examples 
with  deterministic  and  stochastic  permeability  are  given  in  Section  5.  Finally, 
concluding  remarks  are  provided  in  Section  6. 


2  Problem  definition 


In  this  section,  we  follow  the  notation  in  [42] .  Let  us  define  a  complete  probabil¬ 
ity  space  (D,  JF,  V)  with  sample  space  D  which  corresponds  to  the  outcomes  of 
some  experiments,  JF  C  2^  is  the  cr- algebra  of  subsets  in  fl  and  P  :  JF  ^  [0, 1] 
is  the  probability  measure.  Also,  let  us  define  D  as  a  d-dimensional  bounded 
domain  D  C  {d  =  2,3)  with  boundary  dD.  The  governing  equations  for 
immiscible  and  incompressible  two-phase  flow  in  porous  media  consists  of  an 
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elliptic  equation  for  fluid  pressure  and  a  transport  equation  for  the  movement 
of  fluid  phases.  For  simplicity,  we  will  neglect  the  effects  from  gravity,  capil¬ 
lary  forces  and  assume  that  the  porosity  is  a  constant.  The  two  phases  will  be 
referred  to  as  water  and  oil,  denoted  as  w  and  o,  respectively.  The  total  Darcy 
velocity  u  and  the  pressure  p  satisfy  for  P-almost  everywhere  (a.e.)  in  Vt  the 
following  SPDEs  [17] 

V  ■  u  =  q,  u  = —K{x,uj)Xt\/p,  Vcc  G  D,  (1) 

with  the  following  boundary  conditions 

p  =  p  on  dDp,  u  n  =  u  on  dDu-  (2) 

The  total  velocity  w  =  Wo  +  is  a  sum  of  the  velocities  of  oil  Uo  and  water 
Uw  q  is  a  volumetric  source  term  which  is  assumed  0  throughout  the  paper. 
The  random  permeability  tensor  K  is  assumed  to  be  diagonal  and  uniformly 
positive  definite.  In  addition,  we  will  assume  K  is  a  stochastic  scalar  function. 
The  total  mobility  is  given  by  Aj  =  -|-  Ao,  where  A,  models  the  reduced 

mobility  of  phase  i  due  to  the  presence  of  the  other  phase.  Without  loss  of 
generality,  we  assume  that  the  boundary  conditions  are  deterministic  and  that 
the  Neumman  condition  is  homogeneous,  w  =  0  on  dDu. 

Furthermore,  to  assess  the  quality  of  the  multiscale  model,  the  unit  mobility 
ratio  displacement  model  is  used,  i.e.  A^  =  5”,  Aq  =  1  —  5”  and  hence  Aj  =  1, 
where  S  is  the  water  saturation.  Under  these  assumptions,  the  water  saturation 
equation  is  given  by 

dS{x  ^,  ^)  +  .  \/S{x,  t,  Cj)  =  0,  yxeD,te[0,T].  (3) 


Since  the  permeability  77  is  a  stochastic  function,  all  the  unknowns  p,  u  and  S 
are  also  stochastic.  Therefore,  our  complete  stochastic  model  is:  Find  stochas¬ 
tic  functions  and  S'  :  D  x  [0,  T]  x  D  ^  E  for 

"P-almost  every  where  (a.e.)  a;  G  D  such  that  the  following  equations  hold: 


V  •  tr(£c,  a;)  =  0,  u{x,u})  = —K{x,u})Vp{x,u})  \/x  E  D,  (4) 

^^^^^^^  +  u{x,t,u)-VS{x,t,Lj)  =0,  yxeD,te[0,T],  (5) 

at 

with  the  boundary  conditions 

p  =  p  on  dDp,  u  n  =  0  on  dDu,  (6) 

together  with  appropriate  initial  and  boundary  conditions  for  S.  Computa¬ 
tion  with  this  model  is  much  more  efficient  than  using  the  actual  two-phase 
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flow  model  because  the  pressure  and  saturation  equations  are  effectively  de¬ 
coupled.  Throughout  this  paper,  the  Darcy  velocity  u  is  first  computed  using 
the  mixed  finite  element  heterogeneous  multiscale  method  developed  in  Sec¬ 
tion  4.1  and  then  the  saturation  equation  is  solved  using  a  upwinding  finite 
element  scheme  [51]  in  Section  4.2.  Although  these  equations  differ  from  the 
actual  flow  equations,  they  do  capture  many  important  aspects  of  two-phase 
flow  problems.  Specifically,  the  effects  of  the  heterogeneity  are  often  similar  in 
the  unit  mobility  and  two- phase  flow  problems  [52]. 


2.1  The  finite- dimensional  noise  assumption  and  the  Karhunen-Loeve  ex¬ 
pansion 


We  employ  the  ‘finite-dimensional  noise  assumption’  [39]  and  using  the  Karhunen- 
Loeve  (K-L)  expansion  [53]  we  approximate  any  second-order  stochastic  pro¬ 
cess  with  a  finite-dimensional  representation. 

Geostatistical  models  often  suggest  that  the  permeability  field  is  a  weakly  or 
second-order  stationary  random  field  such  that  the  mean  log-permeability  is 
constant  and  its  covariance  function  only  depends  on  the  relative  distance  of 
two  points  rather  than  their  actual  location  [7].  Denote  =  log(K)  and 

its  covariance  function  by  Rg(xi,  X2),  where  Xi  and  X2  are  spatial  coordi¬ 
nates.  By  definition,  the  covariance  function  is  real,  symmetric,  and  positive 
definite.  All  its  eigenfunctions  are  mutually  orthonormal  and  form  a  complete 
set  spanning  the  function  space  to  which  G{x,  cu)  belongs.  Then  the  truncated 
K-L  expansion  takes  the  following  form: 


G{x,  cj)  =  E[G'(a;)]  + 

i=l 


(7) 


where  are  uncorrelated  random  variables.  Also,  fiix)  and  Aj  are 

the  eigenfunctions  and  eigenvalues  of  the  correlation  function,  respectively. 
They  are  the  solutions  of  the  following  eigenvalue  problem: 

/  RG{Xi,X2)(t)i{x2)dX2  =  Xifiixi).  (8) 

JD 

The  number  of  terms  needed  to  approximate  a  stochastic  process  depends  on 
the  decay  rate  of  the  eigenvalues.  Generally,  a  higher-correlation  length  would 
lead  to  a  rapid  decay  of  the  eigenvalues. 


When  using  the  K-L  expansion,  we  here  assume  that  we  obtain  a  set  of  mu¬ 
tually  independent  random  variables.  Denote  the  probability  density  func¬ 
tions  of  as  Pi,  i  =  Let  T,  be  the  image  of  1^.  Then 

p(Y)  =  Hill  Pi{Yi)  is  the  joint  probability  density  of  Y  =  (Yi,  •  •  •  ,  Yjv)  with 
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support  r  =  rixr2X---xrjvG  IR'^.  Then  the  stochastic  log  permeability 
can  be  represented  by  G{x,u})  =  G{x,Yi, . . .  ,Yn)  =  G{x,Y). 


2.2  Stochastic  variational  formulation 


By  using  the  Doob-Dynkin  lemma  [42],  the  solutions  of  Eqs.  (4)  and  (5)  can 
be  described  by  the  same  set  of  random  variables  Following  [48], 

we  define  appropriate  function  spaces  that  encode  variations  of  the  function 
in  the  physical  domain  V  and  in  the  stochastic  space  F. 

In  the  physical  space,  we  introduce  the  following  common  functional  spaces  [16,48] 
W  =  L‘^{D)  =  i^p:  \p\‘^dx  =11  p  ||i2p)<  +oo| ,  (9) 

with  inner  product 

{p,q)  =  {p,q)L^D)  ■=  j^P  q  dx,  p,qeL‘^{D),  (10) 

and 

H{div,  D)  =  {u:ue  Vue  ,  (11) 

with  inner  product 

{u,v)  =  {u,v)H(div,D) J^U  ■  V  dx,  u,v  e  H{dw,D).  (12) 
We  will  also  make  use  of  the  following  space: 

V  =  iFo,u(div,  D)  =  {u  :  u  e  H{div,  D),  u  ■  n  =  0}  .  (13) 

The  duality  product  is  defined  as: 

{u,p)  =  {u,p)oDp  ■■=  [  updx,  ueH^^'^{D),  peH~^/^{D).  (14) 

JdDp 

The  functional  space  in  T  is  defined  as  follows: 

u  =  Ll(T)  =  |p :  [f  |p(r)|V(y)<iy)  c»| .  (is) 

By  taking  its  tensor  product  with  the  previous  deterministic  spaces,  one  can 
form  the  stochastic  functional  spaces: 

yp  =  ZY0lT,  V  =  U®V.  (16) 

Multiplication  of  Eqs.  (4)  and  (5)  by  appropriate  test  functions  and  integration 
by  parts  leads  to  the  following  weak  formulations:  Find  u  e  7i,  p  e  W  such 
that 


v)p{Y)dY  -  ^(V  •  v,p)p{Y)dY 

=  -  J^{vn,p)p{Y)dY,  yvev, 

[  {i,\/  ■u)p{Y)dY  =  o,  View, 


and  S'  e  W  for  each  t  E  [0,  T]  such  that 


X 


,  ^ .  n  p(y )<iy  +  y^(u .  vs,  q)p{Y)dY = o.  v«  e  w. 


(17) 

(18) 


(19) 


We  assume  without  loss  of  generality  that  the  support  of  the  random  variables 
yj  is  P  =  [0, 1]  for  z  =  1,  •  •  •  ,  and  thus  the  bounded  stochastic  space  is  a 
A/’-hypercube  F  =  [0, 1]'^,  since  any  bounded  stochastic  space  can  always  be 
mapped  to  the  above  hypercube. 


3  Adaptive  sparse  grid  collocation  method  (ASGC)  and  High  di¬ 
mensional  model  representation  technique  (HDMR)  for  the  so¬ 
lution  of  SPDEs 


The  original  infinite-dimensional  stochastic  problem  is  now  restated  as  a  fi¬ 
nite  A-dimensional  problem.  Then  we  can  apply  any  stochastic  method  in 
the  random  space  and  the  resulting  equations  become  a  set  of  deterministic 
equations  in  the  physical  space  that  can  be  solved  by  any  standard  deter¬ 
ministic  discretization  technique,  e.g.  the  finite  element  method.  The  solution 
to  the  above  SPDEs  Eqs.  (17)-(19)  can  be  regarded  as  stochastic  functions 
taking  real  values  in  the  stochastic  space  F.  For  example,  we  can  consider  the 
pressure  as  a  stochastic  function  p  :  F  — >  M  and  we  use  the  notation  p{Y) 
to  highlight  the  dependence  on  the  randomness.  Then  it  can  be  shown  that 
the  weak  formulation  Eqs.  (17)-(19)  is  equivalent  to  [38]:  for  a.e.  p  G  F  the 
following  deterministic  weak  form  equations  hold: 


{K  ^u,v)  —  {p,^/ ■  v)  = —{p,v  •  n),  \/veV 


{i,v-u)  =  o,  yiew 


'dS 


—  ,q\+{q,u-VS)  =  0,,  MqeW 


(20) 

(21) 

(22) 


This  nature  is  utilized  by  the  stochastic  collocation  method  which  constructs 
the  interpolant  of  the  stochastic  function  in  F  using  only  the  solutions  to  the 
above  deterministic  problems  at  chosen  sample  points. 
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3.1  Adaptive  sparse  grid  colloeation  method  (ASGC) 


In  this  section,  we  briefly  review  the  development  of  the  ASGC  strategy.  For 
more  details,  the  interested  reader  is  referred  to  [42], 

The  basic  idea  of  this  method  is  to  employ  a  flnite  element  approximation  for 
the  spatial  domain  and  approximate  the  multi-dimensional  stochastic  space 
r  using  interpolating  functions  on  a  set  of  collocation  points  G  F. 

Suppose  we  can  And  a  flnite  element  approximate  solution  S  to  the  deter¬ 
ministic  solution  of  the  problem  in  Eqs.  (20)-(22),  we  are  then  interested  in 
constructing  an  interpolant  of  S  by  using  linear  combinations  of  the  solutions 
The  interpolation  is  constructed  by  using  the  so  called  sparse  grid 
interpolation  method  based  on  the  Smolyak  algorithm  [42],  In  the  context  of 
incorporating  adaptivity,  we  have  chosen  the  collocation  points  based  on  the 
Newton-Cotes  formulae  using  equidistant  support  nodes.  The  corresponding 
basis  function  is  the  multi-linear  basis  function  constructed  from  the  tensor 
product  of  the  corresponding  one-dimensional  functions. 

Any  function  /  :  D  x  F  ^  M  can  now  be  approximated  by  the  following 
reduced  form: 

!(x,Y)=  Y.  EdW-“j(V),  (23) 

||i||^A(+r  j 

where  the  multi- index  i  =  {ii, ...  An)  G  the  multi-index  j  =  (ii, . . . ,  jA)  G 
and  ||i||  =  ii  +  ■  ■  ■  -\-  in-  r  is  the  sparse  grid  interpolation  level  and  the 
summation  is  over  collocation  points  selected  in  a  hierarchical  framework  [42] . 
Here,  tcj  is  the  hierarchical  surplus,  which  is  just  the  difference  between  the 
function  value  at  the  current  point  and  interpolation  value  from  the  coarser 
grid.  The  hierarchical  surplus  is  a  natural  candidate  for  error  control  and 
implementation  of  adaptivity. 

After  obtaining  the  expression  in  Eq.  (23),  it  is  also  easy  to  extract  the  statis¬ 
tics  [42].  The  mean  of  the  random  solution  can  be  evaluated  as  follows: 

E[/(a^)]=  E  [  a\{Y)dY,  (24) 

||i||^Ar+r  j 

where  the  probability  density  function  p{Y)  is  1  since  the  stochastic  space 
is  a  unit  hypercube  [0, 1]^.  As  shown  in  [42],  the  multi-dimensional  integral 
is  simply  the  product  of  the  ID  integrals  which  can  be  computed  analyti¬ 
cally.  Denoting  J  a^^{Y)dY  =  /j,  we  can  rewrite  Eq.  (24)  as  Eg  [/ (a?)]  = 

E||i||<v+rEjWj(*)  •  ly 
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We  now  define  the  error  indicator  as  follows: 


7i  = 


w\{x)  •  /] 


L2(D) 


E||i||-iv-i[/] 


L2(D) 


(25) 


Here,  the  L2  norm  is  defined  in  the  spatial  domain.  This  error  indicator  mea¬ 
sures  the  contribution  of  each  term  in  Eq.  (24)  to  the  integration  value  (mean) 
relative  to  the  overall  integration  value  computed  from  the  previous  interpo¬ 
lation  level.  In  addition  to  the  surpluses,  it  also  incorporates  information  from 
the  basis  functions.  This  makes  the  error  7]  to  decrease  to  a  sufficient  small 
value  for  a  large  interpolation  level.  Therefore,  for  a  reasonable  error  thresh¬ 
old,  this  error  indicator  guarantees  that  the  refinement  would  stop  at  a  certain 
interpolation  level. 


The  basic  idea  of  adaptive  sparse  grid  collocation  (ASGC)  method  here  is  to 
use  the  error  indicator  7]  to  detect  the  smoothness  of  the  solution  and  refine 
the  hierarchical  basis  functions  aj  whose  magnitude  satisfies  7]  ^  where  £ 
is  a  predefined  adaptive  refinement  threshold.  If  this  criterion  is  satisfied,  we 
simply  add  the  2N  neighbor  points  of  the  current  point  to  the  sparse  grid  [42] . 

However,  this  method  still  suffers  from  the  curse  of  dimensionality.  In  the  next 
section,  a  novel  dimension  decomposition  method  is  reviewed  to  transform  the 
A-dimensional  problem  into  several  low-dimensional  sub-problems. 


3.2  High  dimensional  model  representation  (HDMR) 


In  this  section,  the  basic  concepts  of  HDMR  are  briefly  reviewed  following 
closely  the  notation  in  [43] .  For  a  detailed  description  of  the  theory  applied  to 
stochastic  systems,  the  interesting  reader  may  refer  to  [43]. 

In  order  to  alleviate  the  curse  of  dimensionality,  we  have  combined  ASGC 
with  the  adaptive  stochastic  high  dimensional  model  representation  (HDMR) 
technique.  HDMR  represents  the  model  outputs  as  a  finite  hierarchical  corre¬ 
lated  function  expansion  in  terms  of  the  stochastic  inputs  starting  from  lower- 
order  to  higher-order  component  functions.  In  that  work,  the  CUT-HDMR  is 
adopted  to  construct  the  response  surface  of  the  stochastic  solution.  Within 
the  framework  of  CUT-HDMR,  a  reference  point  Y  =  ^Ui,  F2,  •  •  • ,  Ujv)  is 
first  chosen.  According  to  our  past  experience,  the  mean  vector  of  the  random 
input  W  is  a  good  choice  for  the  reference  point.  Then  HDMR  is  given  in  a 
compact  form  as  [43] 

/(V)  =  E  E(-l)'“''''''/(Vv)lr.F\y'..  (26) 

uCV  vCu 
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for  a  given  set  u  C  I),  where  V  :=  A^}  denotes  the  set  of  coordinate 

indices  and  we  define  fiY^)  =  fiY)-  Here,  l^v  denotes  the  |v|-dimensional 
vector  containing  those  components  of  Y  whose  indices  belong  to  the  set  v, 
where  |v|  is  the  cardinality  of  the  corresponding  set  v,  i.e.  l^v  =  (b^)iev-  The 
notation  Y  =  Y  \  Yy,  means  that  the  components  of  Y  other  than  those 
indices  that  belong  to  the  set  v  are  set  equal  to  those  of  the  reference  point. 
For  example,  if  v  =  {1,3,5},  then  |v|  =  3  and  /(l^v)  is  a  function  of  only 
three  random  variables  ^1,13,15  while  the  other  dimensions  satisfy  =  Fj 
for  i  e  P  and  i  ^  v. 

Therefore,  the  A^-dimensional  stochastic  problem  is  transformed  to  several 
lower-order  |v|-dimensional  problems  /(Fv)  which  can  be  easily  solved  by 
the  ASGC  as  introduced  in  the  last  section: 

/(y)=  ^  5;<J(x).oj(rv),  (27) 

uCVyCu  ||i||^A^+r  j 

where  ||i||  =  ii  +  ■  ■  ■  -\-  i|v|,  w^{x)  are  the  hierarchical  surpluses  for  different 
sub- problems  indexed  by  v  and  aj(Fv)  is  only  a  function  of  the  coordinates 
which  belong  to  the  set  v.  It  is  noted  that  the  interpolation  level  r  may  be 
different  for  each  sub-problem  according  to  their  regularity  along  the  particular 
dimensions  which  is  controlled  by  the  error  threshold  e. 

In  addition,  it  is  also  easy  to  extract  statistics  as  introduced  in  Section  3.1  by 
integrating  directly  the  interpolating  basis  functions.  Let  us  denote 

(28) 

vCu  llill^jv+r  j 

as  the  mean  of  the  component  function  /„.  Then  the  mean  of  the  HDMR 
expansion  is  simply  E  [/  (F)]  =  Ju-  To  obtain  the  variance  of  the  solu¬ 

tion,  we  can  similarly  construct  an  approximation  for  v?  and  use  the  formula 
Var  [u  (x)]  =  E  [u^(a?)]  —  (E  [u  (a?)])^. 

It  is  noted  that  the  solution  method  of  each  sub-problem  is  not  limited  to 
ASGC.  It  is  also  possible  to  use  the  sparse  grid  based  on  Gauss  quadrature 
rule  to  integrate  the  component  functions  of  the  CUT-HDMR  in  order  to 
obtain  the  mean  and  the  standard  deviation  directly.  In  this  case,  Eq.  (28) 
can  be  rewritten  as 

j„=  i:(-i)i“Hvi  ^  5;;ny/(x,yj)  (29) 

vCu  ||i||<JV-|-r  j 

where  is  the  quadrature  weight  and  f{x,  F)  is  the  function  value  at  the  col¬ 
location  points.  The  advantage  of  this  method  is  its  higher  accuracy  than  the 
linear  interpolation.  However,  it  does  not  provide  the  function  approximation. 
In  this  work,  we  would  also  like  to  extend  our  previous  ASGC  formulation  to 
include  the  Gauss  quadrature  rule. 
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As  shown  in  [43],  it  is  not  necessary  to  compute  all  the  terms  in  Eq.  (26).  We 
can  take  only  a  subset  S  of  all  indices  u  CD  while  retaining  the  approxima¬ 
tion  accuracy.  Therefore,  we  can  define  an  interpolation  formula  Asf  for  the 
approximation  of  /  which  is  given  by 


^5/:=E^(/u)-  (30) 

Here,  A{fu)  is  the  sparse  grid  interpolant  of  the  component  function  /u  and 
Asf  is  the  interpolant  of  the  function  /  using  the  proposed  method  with  the 
index  set  S.  It  is  common  to  refer  to  the  terms  {/u  :  |u|  =  m}  collectively  as 
the  “order-m  terms” .  Then  the  expansion  order  for  the  decomposition  Eq.  (30) 
is  defined  as  the  maximum  of  m.  Note  that  the  number  of  collocation  points 
in  this  expansion  is  defined  as  the  sum  of  the  number  of  points  for  each  sub¬ 
problem  from  Eq.  (27),  i.e.  M  = 


Based  on  this,  we  have  also  developed  the  adaptive  version  of  HDMR  to  find 
the  optimal  set  <S.  A  weight  is  defined  for  each  expansion  term  as: 


Vu  = 


\\Ju\\L2iD) 


Sve5,|v|<|u|-1 


L2(D) 


(31) 


We  always  construct  the  zeroth-  and  first-order  HDMR  expansion  where  the 
computational  cost  is  affordable  even  for  very  high-dimensions  and  only  select 
those  terms  whose  weight  is  greater  than  a  predefined  error  threshold  6i .  Then 
we  define  the  important  dimensions  as  those  whose  weights  are  larger  than  a 
predefined  error  threshold  9i.  Now,  the  set  T>  in  Eq.  (26)  only  contains  these 
important  dimensions  instead  of  all  the  dimensions.  However,  not  all  the  pos¬ 
sible  terms  are  computed.  Instead,  we  adaptively  construct  higher-order  com¬ 
ponent  functions  increasingly  from  lower-order  to  higher-order.  Similarly,  the 
important  component  functions  are  defined  as  those  whose  weights  are  larger 
than  the  predefined  error  threshold  6i.  We  put  all  the  important  dimensions 
and  higher-order  terms  into  a  set  T,  which  is  called  the  important  set.  When 
adaptively  constructing  HDMR  for  each  new  order,  we  only  calculate  the  term 
/u  whose  indices  satisfy  the  admissibility  relation: 


u  E  V  and  v  C  u  ^  v  G  T. 


(32) 


In  other  words,  among  all  the  possible  indices,  we  only  want  to  find  the  terms 
which  can  be  computed  using  the  previous  known  important  component  func¬ 
tions.  In  this  way,  we  find  those  terms  which  may  have  significant  contribution 
to  the  overall  expansion  while  ignoring  other  trivial  terms  thus  reducing  the 
computational  cost  for  high-dimensional  problems.  In  addition,  if  the  relative 
error  between  two  consecutive  orders  is  smaller  than  another  threshold  62, 
the  HDMR  expansion  is  considered  converged  and  the  construction  stops.  Eor 
more  details,  please  refer  to  [43]. 
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4  Spatial  finite  element  discretization 


As  stated  in  Section  3,  in  order  to  ntilize  the  HDMR  Eq.  (27),  we  only  need 
to  seek  the  solution  (m,p,  S)  at  each  collocation  point  in  the  stochastic  space 
r.  In  other  words,  our  goal  is  reduced  to:  for  each  permeability  realization 
R'W(a;)  =  K{x,Yi),i  =  we  solve  the  deterministic  problem:  find 

y^b)  £  pb)  £  \Y  and  G  W  such  that  for  i  =  1, . . . ,  M 


\v)-{p^\V-v)  =  - 

{p,v  ■  n),  Vv  e  V, 

(33) 

(/,V-w®)  =  0, 

Vt  G  IT, 

(34) 

at 

+  ■  VA(*))  =  0, 

Vg  G  IT. 

(35) 

In  this  section,  mixed  finite  element  methods  are  introduced  to  solve  the  above 
equations  in  the  spatial  domain.  Since  the  pressure  Eqs.  (33)  and  (34)  are 
effectively  decoupled  from  the  saturation  Eq.  (35),  we  will  first  introduce  the 
multiscale  method  to  find  u,p  and  then  use  the  upwinding  finite  element 
method  to  find  S.  To  simplify  the  notation,  we  will  omit  the  superscript  i  and 
assume  the  deterministic  equations  are  satisfied  for  an  arbitrary  permeability 
sample  in  the  stochastic  space. 


4-1  Mixed  finite  element  heterogeneous  multiseale  method  (MxHMM) 


In  the  porous  media  flow  problem,  the  heterogeneity  of  the  permeability  field 
will  have  a  great  impact  on  the  global  flow  conditions.  In  order  to  resolve 
the  fine-scale  velocity  accurately  with  lower  computational  cost,  a  multiscale 
method  is  needed.  In  addition,  the  mixed  finite  element  method  is  also  re¬ 
quired  to  compute  the  velocity  and  pressure  simultaneously,  if  we  want  to 
have  an  accurate  velocity  and  ensure  mass  conservation.  We  can  identify  at 
least  two  main  multiscale  methods:  multiscale  finite  element  or  finite  volume 
methods  [7,9,11]  and  the  variational  multiscale  methods  [14,16].  In  this  sec¬ 
tion,  we  will  develop  a  new  multiscale  method  which  is  based  on  the  framework 
of  the  heterogeneous  multiscale  method  [22] .  We  present  the  discretization  and 
methodology  for  a  two-dimensional  system.  Extension  to  three-dimensions  is 
straightforward . 

Consider  a  partition,  Th  fo  the  domain  D  into  non-overlapping  elements  ej, 
Th  =  where  Nh  is  the  number  of  elements  of  the  grid.  Define  also 

the  skeleton  of  the  partition,  SVh  =  U^i  where  Mh  is  the  number  of 
element  faces  denoted  by  Pa-  The  partition  Th  is  regarded  as  the  fine-scale  grid. 
The  multiscale  permeability  is  defined  as  a  cell- wise  constant  on  this  grid.  To 
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implement  the  multiscale  method,  we  also  consider  a  coarse-scale  partition 
of  the  same  domain  D.  Denote  this  partition  as  %  =  \J^i  Ei.  Denote  by 
SVc  =  U^i  Aa  the  associated  skeleton  of  the  coarse-scale  discretization.  Here, 
Nc  is  the  number  of  coarse  elements  and  M^.  is  the  number  of  coarse  element 
faces  denoted  by  A^.  In  order  to  conserve  the  mass  at  the  coarse-scale,  we  also 
assume  for  simplicity  that  the  partitions  7^  and  7^  are  nested,  conforming 
and  consist  of  rectangular  elements.  Fig.  1  shows  a  hne  grid  (finer  lines)  and 
a  corresponding  coarse  grid  (heavier  lines). 


Fig.  1.  Schematic  of  the  domain  partition:  (a)  fine-  and  coarse-scale  grids,  (b) 
fine-scale  local  region  in  one  coarse  element. 

Now  consider  the  finite  dimensional  subspaces  on  the  coarse-scale  Vc  E  V  and 
Wc  G  W.  The  mixed  finite  element  method  approximation  of  Eqs.  (33)- (34) 
on  the  coarse-scale  reads:  Find  the  coarse-scale  (ucPc)  G  14  x  H4  such  that 


(K  ^Uc,Vc)  -  (Pc,'^  ■Vc)  =  -{po,Vc-n),  V  •Uc  G  14,  (36) 

(4,V-We)  =  0,  V/eGbFe.  (37) 

Note  that  14  and  H4  should  satisfy  the  discrete  inf-sup  condition  [54].  In  this 
work,  14  is  taken  to  be  the  lowest-order  Raviart- Thomas  space  [8],  RTq{Tc) 
and  H4  is  taken  to  be  the  space  of  piece-wise  constants  on  the  coarse-scale 
mesh,  Pq{%.)-  Other  choices  can  be  found  in  [54].  Therefore,  we  define  the 
finite  element  space  for  the  coarse-scale  velocity  as: 

f  'I 

I4=|wc:'Ue  =  ^7/.X,  <  =  0,  VA„GaD„|,  (38) 

where  is  the  RTq  basis  functions  on  the  uniform  mesh  of  rectangular  el¬ 
ements  associated  with  the  coarse  element  face  A^.  For  a  reference  element 
E  =  [x{,x^]  X  [x2,X2]  with  the  area  \E\,  there  are  four  vector  RTq  basis 
functions  with  non-zero  support: 
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^1= 

{xf  —  xi)/{xf  —  Xi  ),  0 

T 

,  ^2  = 

0,  (xf  -  X2)/(xf  -  X2) 

^3  = 

[(xi  -xf)/(xf -xf),ol 

T  1 

>  '04  = 

0,  (X2  -  x\)l{x^  -  X2  )] 

5 


(39) 

(40) 


The  basis  functions  satisfy  the  properties  such  that  -0^  •  Uj  =  1  ii  i  =  j, 
otherwise  •  Uj  =  0  for  hi  =  1, . . . , 4.  Therefore  is  value  of  the  coarse- 
scale  flux  at  the  middle  point  of  the  side  A^,  i.e.  Uc  -  ria  =  where  is  the 
unit  outer  normal  to  the  interface  A^.  The  coarse-scale  pressure  approximation 
is  piecewise  constant  on  the  coarse- mesh  and  Po{Tc)  is 


Wc=  Pc 


(41) 


where  is  the  coarse-scale  pressure  basis  function  for  the  coarse  element  i 
defined  as 

1,  if  a;  e  Ei, 

0,  if  a;  ^  Ei. 

Pi  is  the  corresponding  pressure  degree  of  freedom  (the  average  pressure  in 
coarse  element  Ei). 


It  is  obvious  that  all  the  fine-scale  information  is  included  in  the  bilinear  form 
{K~^Uc^Vc)-  Denote  A  —  {Aij)  the  global  matrix  for  the  bilinear  form,  where 

Ai  =  'fPiix)  •  (43) 

We  could  evaluate  Eq.  (43)  by  the  2x2  Gauss  quadrature  rule:  let 

fij{x)  ^il)^i{x)  ■  K-\x)il)‘j{x),  (44) 

then 


where  and  Tk,k  =  1, . . . ,  4  are  the  quadrature  points  and  weights  (including 
the  determinant  of  the  Jacobian  matrix)  in  the  coarse  element  E,  respectively. 
It  is  obvious  that  any  realization  of  the  permeability  field  at  the  quadrature 
point  K  (^^)  is  not  able  to  capture  the  full  information  at  the  subgrid  scale  in 
the  coarse  element  since  the  size  of  the  coarse  element  is  much  larger  than  the 
characteristic  length  scale  of  the  multiscale  permeability  field.  Therefore,  we 
need  to  modify  the  bilinear  form  Eq.  (44)  at  the  quadrature  point  following 
the  framework  of  the  heterogeneous  multiscale  method  [21,26]  as: 

/p(^fc)  =  — br  /  Uik{x)  ■  K~'^Ujk{x)  dx,  A:  =  l,  ...,4,  (46) 

Es^.\ 
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where  Uik{x),  i  =  1, . . . ,  4  is  the  solution  to  the  following  local  subgrid  problem 
in  the  sampling  domain  Es^.  C  fc  =  1, . . . ,  4: 

V  •  iiikix)  =  0,  Uik{x)  =  -KVpikix),  y  X  e  Es^,  (47) 

with  appropriate  boundary  conditions  which  we  will  discuss  below.  Pi(x)  can 
be  considered  as  the  subgrid  pressure. 

First,  we  will  discuss  the  choice  of  the  sampling  domain  Es^  of  the  subgrid 
problem.  In  the  original  problem  definition  of  the  FeHMM  [21,26],  the  coeffi¬ 
cient  of  the  elliptic  equation  (here  K)  is  assumed  to  be  periodic.  Therefore,  the 
sampling  domain  was  taken  around  each  quadrature  point  as  Es^.  =  +  SI, 

where  I  =  (—1/2, 1/2)^  and  S  is  equal  to  one  period  of  the  coefficient  in  the 
elliptic  equation,  as  in  Fig.  2(a).  However,  in  general,  the  permeability  is  not 
periodic.  If  the  sampling  domain  is  too  small,  one  cannot  capture  enough  in¬ 
formation  on  the  subgrid  scale.  According  to  the  numerical  results  in  [55] ,  the 
larger  the  size  of  the  sampling  domain  is,  the  more  accurate  the  computed 
result  is.  Therefore,  we  would  like  to  take  the  sampling  domain  to  be  the  same 
as  the  coarse  element,  i.e.  Es^  =  E,dEs^  =  dE  as  in  Fig.  2(b).  In  addition, 
we  also  assume  that  the  fine  grid  within  each  coarse  element  is  the  same  as 
the  fine-scale  grid  7/,  where  the  permeability  is  defined,  see  Fig.  1(b).  In  this 
way,  we  can  ensure  global  continuity  of  the  flux  across  the  coarse  element. 


* 

•  • 

•  • 

II 

•  • 

•  • 

•  • 

• 

• 

• 

• 

(a)  (b) 

Fig.  2.  (a)  Schematic  of  the  original  HMM  method,  where  the  sampling  domain 
is  around  the  quadrature  point,  (b)  Schematic  of  the  proposed  MxHMM  method, 
where  the  sampling  domain  is  the  same  as  the  coarse  element. 

Remark  1.  Unlike  the  mixed  multiscale  finite  element  method  [9],  where 
the  subgrid  problem  is  limited  to  the  coarse  element,  the  advantage  of  the 
heterogeneous  multiscale  method  here  is  that  the  sampling  domain  is  not 
limited  to  the  coarse  element.  In  fact,  it  can  be  chosen  arbitrarily  to  include 
as  many  coarse  elements  as  necessary.  However,  in  the  present  work  we  still 
solve  the  subgrid  problem  in  only  one  coarse  element.  The  effect  of  the  size  of 
the  sampling  domain  is  reserved  for  later  work. 

Hence,  all  the  subgrid  problems  are  solved  within  the  same  coarse  element. 
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The  only  difference  is  the  applied  boundary  condition.  The  boundary  condi¬ 
tion  of  the  problem  in  Eq.  (47)  plays  a  significant  role  in  the  accuracy  of  the 
multiscale  method  as  discussed  in  [55],  where  three  different  boundary  con¬ 
ditions  are  considered:  the  periodic  boundary  condition,  Dirichilet  boundary 
condition,  and  the  Neumann  boundary  condition.  However,  when  mixed  finite 
element  formulation  is  used  on  the  coarse-scale,  only  the  Neumann  bound¬ 
ary  condition  is  applicable  here.  In  [26],  the  following  Neumann  boundary 
condition  is  proposed: 


Aik  ■  noE  =  ^Uik)  ■  on  SE, 


(48) 


where  denotes  the  value  of  the  t-th  coarse-scale  RTq  finite  element  basis 

function  at  the  quadrature  point  =  1,...,4  and  uqe  denotes  the  unit 
outer  normal  of  the  coarse  element  boundary  dE.  According  to  the  definition 
of  RTq  basis  function  in  Eqs.  (39)- (40),  this  boundary  condition  applies  a 
uniform  flow  with  magnitude  from  one  side  to  the  opposite  side  while 

keeping  no- flow  conditions  on  the  other  two  sides.  The  example  of  '?/’i(^i) 
is  shown  in  Fig.  3.  However,  this  boundary  condition  only  reflects  the  local 
heterogeneity  structure  within  the  current  coarse  element.  It  does  not  contain 
the  flow  condition  across  the  coarse  element  interface  which  is  often  important 
in  guaranteeing  the  continuity  of  flux  on  the  coarse-scale.  Therefore,  we  would 
like  to  propose  a  new  boundary  condition  which  reflects  the  heterogeneous 
structure  across  the  coarse  element  boundary. 

For  a  fine-scale  element  interface  Ua,  denote  the  two  adjacent  fine-scale  ele¬ 
ments  as  Si  and  e^,  i.e.  t'a  =  e*  fl  ^j-  According  to  two-point  flux  approximation 
finite  volume  method,  if  the  element  interface  is  in  the  ^-direction,  the  element 
interface  transmissibility  in  the  ^-dimension  is  defined  by  [56]: 


T  = 

I'a 


2|z/„ 


Axi  Axi  \  ^ 


(49) 


where  \va\  is  the  length  of  the  interface,  Axi  denote  the  length  of  element  Cj  in 
the  x-coordinate  direction,  and  Ki  is  the  permeability  in  element  e*.  Similar 
expression  can  be  defined  in  the  ^-dimension.  The  fine-scale  transmissibility 
of  interface  Ua  reflects  the  flow  condition  across  elements.  Denote  the  total 
applied  flux  along  the  coarse  element  interface  A  due  to  the  value  of  the  f-th 
coarse-scale  basis  functions  at  the  fc-th  quadrature  point  as 


Qik  =  [  i’tiD  ■nds=  |A|i/j)=(^fc)  •  ua.  (50) 

JA 


In  this  work,  we  consider  rectangular  elements  oriented  along  the  coordinate 
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axes.  Hence,  we  modify  the  boundary  condition  Eq.  (48)  to: 


Tv 

Uik  ■  n|A  =  Qik  ■  — - — r,  on  A  c  dE, 


(51) 


where  Qik  is  defined  in  Eq.  (50)  and  is  the  fine-scale  transmissibility  of 
interface  Ua  C  A  as  defined  in  Eq.  (49).  See  for  example  Qn  in  Fig.  3(b).  In 
the  equation  above,  is  the  fine-scale  transmissibility  of  interface  as  defined 
in  Eq.  (50)  for  an  interface  in  the  y-direction.  When  the  interface  is  in  the 
x-direction,  we  change  the  definition  of  accordingly.  Therefore,  the  sum  of 
the  flux  applied  on  the  fine-scale  element  is  equal  to  the  total  flux  applied  on 
the  same  coarse  element  boundary.  We  just  redistribute  the  total  flux  on  the 
coarse-scale  element  boundary  according  to  the  ability  to  transport  the  flow 
at  the  interface  of  each  fine-scale  element.  This  is  clearly  a  better  choice  for 
boundary  condition  since  it  determines  the  flow  conditions  across  the  inter¬ 
block  boundaries. 


<(^,)=[(xW4)/(xf-xf),o' 

u  n  =  0 


u  n  =  0 


e,A/Z?;h,i 

i=5 

QnTjf,T,\v,\ 

i=5 

QnT^/t,T,\v,\ 

i=5 


Fig.  3.  Schematic  of  different  boundary  conditions,  (a)  The  uniform  boundary  con¬ 
dition,  (b)  The  modified  boundary  condition  where  the  flux  is  scaled  according  to 
the  fine-scale  transmissibilities. 


Finally,  our  subgrid  problem  in  a  coarse-element  E  is  defined  as  follows:  For 
each  quadrature  point  =  1, . . . ,  4,  we  seek  the  solution  Uik  to  the  following 
subgrid  problem  for  each  coarse-scale  RT^  basis  function  '0^,  z  =  1, . . . ,  4: 

Vuikix)  =  0,  Uik{x)  =  -KVpikix),  y  X  e  E,  (52) 

with  the  Neumman  boundary  condition  defined  in  Eq.  (51). 
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For  convenience,  we  will  define  the  corresponding  modified  bilinear  form  as: 
for  any  Uc,Vc  G  K 

4 

AhiK-W:V,)  :=  E  E  ^  /  Ukix)  •  K-^Vk{x)  dx,  (53) 

EeTck=l  1-^1“'® 

where  U  and  V  are  defined  through  the  subgrid  problems.  The  assembly  of 
this  bilinear  form  will  be  detailed  in  Section  4.2.  Therefore,  the  MxHMM  ver¬ 
sion  of  Eqs.  (36)-(37)  on  the  coarse-scale  reads:  Find  the  coarse-scale  (Wc,  Pc)  G 
Vc  X  Wc  such  that 


Ah{K  ■  Vc)  =  -{po,Vc- n),  V  G  K,  (54) 

(/e,V-w4  =  0,  yiceWc,  (55) 

with  the  boundary  condition 

Pc  =  P  on  dDp,  Uc  ■  n  =  0  on  dDu-  (56) 

The  major  difference  between  Eqs.  (36)-(37)  and  Eqs.  (54)-(55)  lies  in  the  bi¬ 
linear  form  Ah{-,-),  which  needs  solution  of  the  local  subgrid  problem  Eq.  (52). 
It  is  through  these  subgrid  problems  and  the  mixed  formulation  that  the  ef¬ 
fect  of  the  heterogeneity  on  coarse-scale  solutions  can  be  correctly  captured. 
Unfortunately,  it  is  not  trivial  to  analyze  this  multiscale  method  in  a  general 
case,  but  convergence  results  have  been  obtained  using  the  homogenization 
theory  in  the  case  of  periodic  coefficients  [26]. 


4.1.1  Solution  of  the  subgrid  problems  and  assembly  of  the  bilinear  form 

In  general,  the  subgrid  problem  Eq.  (52)  can  be  solved  through  the  stan¬ 
dard  or  mixed  finite  element  method.  In  the  present  setting,  since  we  are 
only  interested  in  the  velocity,  the  mixed  finite  element  method  is  preferred. 
Let  Eh  =  %i{E)  denote  the  fine  grid  defined  over  one  coarse  element  E.  As 
mentioned  before,  it  coincides  with  the  fine-scale  grid  Th-  The  subgrid-scale 
velocity  functional  spaces  will  be  defined  on  the  fine  grid  Eh  of  each  coarse 
element: 

f  Me  ] 

Fe=  =  e  [ ,  (57) 

where  Me  is  the  number  of  edges  in  E,  and  the  pressure  space  is  defined 
similarly: 

WE=i^p:p=if'lp'l,  f’lePoiEh)^  (58) 

where  Ne  is  the  number  of  elements  in  E.  It  is  noted  that,  as  the  Neumann 
boundary  conditions  in  Eq.  (51)  are  imposed  on  all  boundaries  of  the  coarse 
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element  E,  an  extra  constraint  must  be  added  to  make  the  subgrid  problem 
well  posed.  In  our  implementation,  the  pressure  is  prescribed  to  0  at  one  of 
the  elements  in  Eh. 

The  mixed  finite  element  method  approximation  of  Eq.  (52)  in  coarse  element 
Ei  on  the  subgrid-scale  grid  reads:  Find  the  subgrid-scale  (u,p)  €  14;.  x  II4. 
such  that 


(K  ^u,v)  -  (p,V  ■  v)  =  0,  VveVEi,  (59) 

(l,V-v)  =  0,  VIeWe,,  (60) 

with  the  boundary  condition  Eq.  (51).  It  is  noted  that  for  each  coarse  element, 
we  need  to  solve  4  (number  of  quadrature  points)  x  4  (number  of  basis  func¬ 
tions)  =  16  subgrid  problems.  However,  the  only  difference  between  them  are 
in  the  boundary  conditions.  Therefore,  we  only  need  to  assemble  the  stiffness 
matrix  once  and  solve  the  same  algebraic  problem  with  different  right  hand 
vectors. 

Following  a  standard  assembly  process  for  the  global  matrix  of  the  coarse- 
scale  bilinear  form  Eq.  (53),  we  compute  the  contribution  Ae  to  the  global 
matrix  associated  with  the  coarse  element  E,  where  is  a  4  x  4  matrix. 
Assume  the  solution  of  the  subgrid  problem  at  the  fc-th  Gaussian  point  can 
be  written  as  Uik  =  z  =  1, . . . ,  4.  We  can  write  all  the  solutions  as 

a  4  X  Ne  matrix,  =  {c^j)  where  the  i-th  row  contains  the  subgrid  solution 
corresponding  to  the  Ath  coarse-scale  basis  function  '0^.  Therefore,  the  value 
of  Ae  from  the  k-th  Gauss  point  can  be  denoted  as  Ag  =  {a!%)ij,  where 

=  ^Cii  f  •  tp'^dx  Cjm-  (61) 

^  ^  I TA  I  J  E 

Denoting  the  bilinear  form  matrix  from  the  snbgrid-scale  problem  as  = 
(0>  C  =  Ie  we  can  write: 

=  (62) 

k  =  l  1^1 


Finally,  we  would  like  to  comment  on  the  solution  of  the  linear  systems  result¬ 
ing  from  the  mixed  finite  element  discretization.  The  linear  system  is  indefinite, 
and  it  is  difficult  to  solve  using  a  common  iterative  method.  In  our  implemen¬ 
tation,  we  use  the  Schnr  complement  matrix  to  solve  the  pressure  first  and 
then  solve  the  velocity  [10].  The  linear  system  is  solved  using  preconditioned 
conjugate  gradient  method.  All  the  implementations  are  based  on  the  data 
structure  of  the  numerical  library  PETSc  [57]. 


21 


4-2  Reconstruction  of  the  fine-scale  velocity  and  solution  of  transport  equa¬ 
tion 


So  far  we  have  described  the  development  of  the  mixed  finite  element  heteroge¬ 
neous  multiscale  method  for  the  solution  of  the  coarse-scale  velocity.  However, 
in  order  to  simulate  the  transport  equation  accurately,  we  need  to  reconstruct 
the  fine-scale  velocity  using  the  coarse-scale  velocity  and  the  subgrid  perme¬ 
ability.  It  is  noted  that  the  coarse-scale  velocity  is  not  conservative  at  the 
fine-scale.  In  order  to  obtain  a  mass-conservative  fine-scale  velocity,  we  solve 
Darcy’s  equation  within  each  coarse  element  E  using  Neumann  boundary  con¬ 
dition  given  by  the  coarse-scale  flux  along  the  coarse-element  boundary.  The 
coarse-scale  fiux,  denoted  by  is  directly  given  as  the  solution  of  the  sys¬ 
tem  of  linear  equations  from  the  coarse-scale  discretization.  That  is,  for  each 
E  E  %,  one  solves  the  fine-scale  velocity  Uh  inside  E  by  [17,56] 

V  ■Uh  =  0,  Uh  =  -KVph,  yx  e  E,  (63) 


with  the  boundary  condition  similar  to  the  one  used  in  Eq.  (51): 


Uh  ■  n\A  =  Q"  • 


T 

-J-  IJ. 


^1^6  I  I 


on  A  C  dE, 


(64) 


where  Q'^  is  the  coarse-scale  flux  across  the  coarse  element  interface  A,  and 
is  the  fine- scale  transmissibility  of  interface  t'a  C  A.  Since  mixed  finite 
element  method  to  solve  the  coarse-scale  equations,  the  coarse-scale  flux  is 
obtained  directly.  Similar  to  the  subgrid  problem,  Neumann  boundary  condi¬ 
tion  is  applied  on  all  the  boundaries  of  the  coarse  element.  To  obtain  a  unique 
solution  of  the  above  problem,  the  pressure  is  fixed  to  the  coarse-scale  pres¬ 
sure  Pc  in  the  center  element  of  the  mesh  Eh-  As  indicated  in  [17,56,58],  this 
reconstruction  step  guarantees  the  continuity  of  the  flux  across  the  fine-scale 
elements  between  two  coarse  blocks  and  accounts  for  subgrid  heterogeneity. 
It  also  forces  the  sum  of  the  fine  grid  fluxes  to  be  equal  to  the  corresponding 
coarse-scale  flux.  In  this  way,  the  resulting  fine-scale  velocity  is  conservative 
on  fine-scale  grid  as  well  as  the  coarse-scale  grid. 


For  the  solution  of  the  saturation  equation,  we  use  the  upwinding  finite  element 
method  [7,51],  which  is  equivalent  to  the  standard  upstream  weighted  finite 
volume  method.  We  also  approximate  the  saturation  as  a  piecewise  constant 
in  each  fine-scale  element  e,  Po{Th),  the  same  as  the  pressure  space.  Given 
the  discrete  reconstructed  fine-scale  velocity  field  Uh,  for  a  fine-scale  element 
e  E  Eh-  We  define  the  inflow  boundary  of  the  element  as  9e_,  if  W/i  •  n  <  0  on 
de_  and  similarly  the  outflow  boundary  as  cle+,  if  it/j  •  n  ^  0  on  cle+.  For  any 
piecewise  constant  function  Sh  over  the  mesh  Th,  the  upwinding  value  on  de 
is  defined  as  Sh  and  is  equal  to  the  interior  trace  of  Sh  if  on  and  equal 
to  the  exterior  trace  of  Sh  if  on  de-.  In  addition,  we  also  assume  =  0  on 
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de_r\dD. 


Therefore,  the  weak  formulation  of  the  upwinding  scheme  is  to  find  Sh  G  Wh 
such  that 


Let  At  be  the  time  step  and  denote  by  the  approximation  of  the  water 
saturation  in  fine-scale  element  ej  at  time  t^.  Then  the  discrete  form  of  the 
saturation  Eq.  (65)  is: 

(66) 

|e*| 

Here  |ej|  is  the  area  of  the  element  e,.  fij{S)  =  max{sign(gjj)5'j,  —sign{qij)Sj} 
is  the  upwinding  water  saturation  for  the  interface  =  dci  fl  dcj.  Finally,  the 
fiux  across  the  boundary  is  qij  =  Uh  •  riij  ds  where  riij  is  the  unit  normal  to 
Vij  pointing  from  e,  to  ej.  It  is  noted  that  in  Eq.  (66),  only  the  fiux  q^j  on  the 
each  interface  is  required.  This  value  is  directly  computed  as  the  solution  from 
our  multiscale  approach.  This  is  the  main  reason  why  the  method  discussed 
here  is  better  than  the  stabilized  conforming  finite  element  method  [59]. 

It  is  emphasized  again  that  we  consider  the  transport  problem  with  unit  mo¬ 
bility  ratio,  so  the  saturation  changes  will  not  affect  the  pressure  or  velocity. 
Therefore,  we  can  first  compute  the  fine-scale  velocity  with  our  multiscale  ap¬ 
proach  and  then  solve  the  transport  equation.  The  fiow  rate  of  produced  oil 
at  the  outlet  boundary  is  denotes  as  qo  and  the  fiow  rate  of  produced  water 
qw  To  assess  the  quality  of  our  multiscale  approach,  we  will  use  the  so  called 
water  cut  curve  F,  which  defines  the  fraction  of  water  in  the  produced  fluid, 
i.e.,  F  =  qw/idw  +  Qo)  as  a  function  of  time  measured  in  pore  volume  injected 
(PVI).  The  water-cut  is  defined  as 

•  'n)S  ds 

where  refers  to  the  part  of  the  boundary  with  outer  fiow,  i.e.  Uh-n  >  0. 

PVI  represents  dimensionless  time  and  is  computed  as 

PVI  =  Jq  dt/Vp,  (68) 

where  Vp  is  the  total  pore  volume  of  the  system,  which  is  equal  to  the  area  of 
the  domain  D  here  and  Q  =  (tt/i  •  n)  ds  is  the  total  fiow  rate. 

The  complete  schematic  of  the  stochastic  multiscale  method  for  porous  media 
fiow  is  illustrated  in  Fig.  4. 


(67) 
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Fig.  4.  Schematic  of  the  developed  stochastic  multsicale  method  for  porous  media 
flow. 

5  Numerical  Examples 


In  the  first  two  examples,  we  solve  the  problem  with  deterministic  permeabil¬ 
ity  in  order  to  validate  the  newly  developed  multiscale  method.  In  the  third 
example,  the  complete  stochastic  problem  with  a  known  covariance  function 
is  addressed. 


5.1  Simulation  in  realistic  two-dimensional  reservoirs 


This  test  case  is  a  two-dimensional  problem  with  a  highly  heterogeneous  per¬ 
meability.  The  permeability  field  shown  in  Fig.  5  is  taken  from  the  top  layer 
of  the  10-th  SPE  comparative  solution  project  [60].  The  fine  grid  on  which 
the  permeability  is  defined  consists  of  60  x  220  gridblocks.  It  has  Dirichlet 
boundary  conditions  p  =  100  on  {x2  =  0},  p  =  0  on  {x2  =  220}  and  Neumann 
boundary  conditions  tt  •  n  =  0  on  both  {xi  =  0},{a:i  =  60}.  We  also  im¬ 
pose  zero  initial  condition  for  saturation  5'(x,  0)  =  0  and  boundary  condition 
S{x,t)  =  1  on  {x2  =  0}. 

The  reference  solution  is  computed  on  the  fine-scale  grid  using  single-scale 
mixed  finite  element  method  directly,  as  shown  in  Fig.  6(a)  and  Fig.  7(a).  We 
also  show  the  solutions  obtained  with  the  MxHMM  method  on  various  coarse 
grids  in  Figs.  6  and  7.  It  is  seen  that  the  fiow  focuses  along  the  region  with 
higher  permeability  while  bypassing  the  low-permeability  areas.  At  the  same 
time,  the  velocity  field  displays  significant  small-scale  structure  correspond- 
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Fig.  5.  Logarithm  of  the  permeability  field  from  the  the  top  layer  of  the  10-th  SPE 
model,  which  is  defined  on  60  x  220  fine  grid. 

ing  to  the  spatial  permeability  variations.  The  multiscale  solution  successfully 
captures  all  the  main  characters  of  the  fine-scale  results  and  compares  very 
well  with  the  fine-scale  solution,  with  the  two  results  being  quite  difficult  to 
distinguish  visually.  As  a  direct  measure  of  the  error  in  the  computed  velocity 


u  ■  u  ,  where  the  corre- 


field,  we  consider  the  norm:  ||t4||2  = 


spending  relative  error  is  given  as  5{u)  =  ||Wref  —  WmsIl/HwrefH-  The  result 
is  given  in  Table  1.  In  general,  the  error  is  larger  with  coarser  grid  which  is 
possibly  due  to  some  large  local  error  in  high  permeability  region  where  the 
velocity  changes  quickly. 

However,  for  reservoir  simulation  the  most  crucial  factor  is  the  transport  prop¬ 
erties  of  a  velocity  field.  That  is,  a  large  local  error  in  the  velocity  field  may 
not  be  crucial  as  long  as  the  overall  transport  properties  are  correct.  There¬ 
fore,  we  give  the  contour  plots  of  the  saturation  at  time  0.4  PVI  for  various 
coarse  grids  in  Fig.  8.  The  four  multiscale  results  compare  very  well  with 
the  reference  solution.  To  assess  the  accuracy  of  the  transport  properties, 
we  measure  the  relative  difference  in  the  saturation  profile  at  a  given  time: 


5{S)  =  (  /  IfSpef— AmsP  dxY^'^  j (  [  |<S'refP  dx)^^'^ .  The  result  is  given  in  Table  1. 
Jd  Jd 


It  is  seen  that  although  the  corresponding  velocity  error  is  larger  for  the  same 
coarse  grid,  the  saturation  error  is  significantly  smaller. 

Finally,  we  consider  the  water  cut,  which  is  shown  in  Fig.  9.  Once  again, 
the  results  compare  well  with  the  reference  solution.  Here,  we  measure  the 
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Fig.  6.  Contour  plots  of  the  x-velocity  component  for  various  meshes:  (a)  60  x  220 
hue-scale  grid,  (b)  30  x  110  coarse  grid,  (c)  15  x  55  coarse  grid,  (d)  10  x  44  coarse 
grid,  (e)  6  x  22  coarse  grid. 


(a)FineScale  (6)30x110  (c)15x55  (ii')10x44  (e)6x22 


Fig.  7.  Contour  plots  of  the  y-velocity  component  for  various  meshes:  (a)  60  x  220 
fine-scale  grid,  (b)  30  x  110  coarse  grid,  (c)  15  x  55  coarse  grid,  (d)  10  x  44  coarse 
grid,  (e)  6  x  22  coarse  grid. 


maximum  error  as  S(F)  =  maxt^ol-fref(^)  “  -^ms(^)|-  The  result  is  shown  in 
Table  1,  where  the  error  is  quite  small.  Note  that  this  is  a  quite  strict  measure, 
since  the  water  cut  curves  tend  to  be  steep  right  after  breakthrough,  and  thus 
a  small  deviation  in  breakthrough  time  may  give  a  large  value  in  the  error 
measure. 

Overall,  through  this  example,  it  is  shown  that  the  introduced  multiscale 
method  is  quite  robust  and  accurate  for  different  mesh  discretizations. 
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(a)  Fine  Scale  (Z?)  30x110  (c)  15x55  (J)  10x44  (e)6x22 


Fig.  8.  Contour  plots  of  Saturation  at  0.4  PVI:  (a)  60  x  220  fine-scale  grid,  (b) 
30  X  110  coarse  grid,  (c)  15  x  55  coarse  grid,  (d)  10  x  44  coarse  grid,  (e)  6  x  22  coarse 
grid. 


PVI 


Fig.  9.  Water  cut  curves  for  various  coarse  grids. 

Table  1 

Relative  errors  for  various  coarse  grids  in  Example  1. 


Errors 

30  X  no 

15  X  55 

10  X  44 

6  X  22 

5{u) 

0.112 

0.159 

0.170 

0.234 

S{S) 

0.025 

0.049 

0.067 

0.124 

5{F) 

0.0033 

0.0019 

0.0101 

0.0165 

5.2  Simulation  in  a  realization  sample  from  a  random  permeability  filed 


In  this  section,  we  consider  only  a  sample  realization  from  a  random  perme¬ 
ability  field,  which  can  be  considered  as  a  deterministic  run  at  a  collocation 
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Fig.  10.  Logarithm  of  the  permeability  held  from  one  sample  of  a  log-normal  per¬ 
meability  hied  dehned  on  100  x  100  hne-scale  grid. 

point  in  a  stochastic  simulation.  The  permeability  is  defined  on  a  100  x  100 
fine-scale  grid,  which  is  shown  in  Fig.  10.  Flow  is  induced  from  left-to-right 
with  Dirichlet  boundary  conditions  p  =  100  on  {xi  =  0},  p  =  0  on  {xi  =  100} 
and  no-flow  homogeneous  Neumann  boundary  conditions  on  the  other  two 
edges.  We  also  impose  zero  initial  condition  for  saturation  5'(x,0)  =  0  and 
boundary  condition  S{x,t)  =  1  on  the  inflow  boundary  {xi  =  0}.  The  refer¬ 
ence  solution  is  again  taken  from  the  single-scale  mixed  finite  element  on  the 
fine-scale  grid  directly.  All  the  errors  are  defined  the  same  as  before. 

In  Figs.  11  and  12,  we  show  the  velocity  contour  plots  of  the  reference  solution 
and  the  multiscale  solution  on  a  25  x  25  coarse  grid.  The  flow  tries  to  go  through 
the  high  permeable  regions  and  bypass  the  low  permeable  regions,  which  is 
clearly  reflected  in  the  saturation  plot  at  time  0.4  PVI  as  shown  in  Fig.  13. 
All  the  three  figures  compare  well  with  the  reference  solutions.  The  relative 
errors  are  shown  in  Table  2.  We  note  the  relatively  small  saturation  errors 
compared  with  the  large  velocity  errors,  which  again  confirms  that  the  large 
local  velocity  errors  may  not  reflect  the  overall  accuracy  of  the  saturation 
results  as  long  as  the  multiscale  method  captures  the  major  feature  of  the 
underlying  permeability  field. 

Water  cut  curves  are  shown  in  Fig.  14  and  the  maximum  error  is  given  in 
Table  2.  All  the  water  cut  curves  are  visually  nearly  the  same.  The  two  de¬ 
terministic  numerical  examples  successfully  validate  the  introduced  multiscale 
model.  Since  the  stochastic  multiscale  framework  only  requires  repeated  solu¬ 
tion  of  the  deterministic  problems  at  different  collocation  points,  it  is  expected 
to  also  have  accurate  statistics  of  the  solution  in  the  stochastic  simulation  as 
shown  in  the  next  example. 
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(a)  Fine  Scale  (b)  25x25 


Fig.  11.  Contour  plots  of  the  x- velocity  component  for  (a)  100  x  100  fine-scale  grid, 
(b)  25  X  25  coarse  grid. 


Fig.  12.  Contour  plots  of  the  ^-velocity  component  for  (a)  100  x  100  fine-scale  grid, 
(b)  25  X  25  coarse  grid. 

Table  2 

Relative  errors  for  various  coarse  grids  in  Example  2. 


Errors 

50  X  50 

25  X  25 

20  X  20 

10  X  10 

5{u) 

0.060 

0.156 

0.183 

0.324 

5{S) 

0.019 

0.065 

0.089 

0.182 

m 

0.0017 

0.0059 

0.0149 

0.0079 

5.3  Simulation  in  random  permeability  field 

In  the  last  two  examples,  we  have  successfully  verified  the  accuracy  of  our 
newly  developed  multiscale  solver.  In  this  example,  we  investigate  the  statis- 
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(a)  Fine  Scale  (^)  25  x  25 


Fig.  13.  Contour  plots  of  Saturation  at  0.4  PVI:  for  (a)  100  x  100  fine-scale  grid, 
(b)  25  X  25  coarse  grid. 


Fig.  14.  Water  cut  curves  for  various  coarse  grids. 

tical  properties  of  the  transport  phenomena  in  random  heterogeneous  porous 
media.  The  domain  of  interest  is  unit  square  [0, 1]^.  Flow  is  still  induced 
from  left-to-right  with  Dirichlet  boundary  conditions  p  =  1  on  {xi  =  0}, 
p  =  0  on  {xi  =  1}  and  no-flow  homogeneous  Neumann  boundary  condi¬ 
tions  on  the  other  two  edges.  We  also  impose  zero  initial  condition  for  satura¬ 
tion  S{x,0)  =  0  and  boundary  condition  S{x,t)  =  1  on  the  inflow  boundary 

{xi  =  0}. 


The  log-permeability  is  taken  as  zero  mean  random  held  with  a  separable 
exponential  covariance  function 

Cov(a;, y)  =  cr^exp 

where  Li  and  L2  are  the  correlation  lengths  in  x  and  y  direction,  respectively. 
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cr  is  the  standard  deviation  of  the  random  field.  The  K-L  expansion  is  used  to 
parameterize  the  field  as 


Y(n;)  =  log  {K{u))  =  ^iUx)Yi, 


i=l 


(70) 


where  the  eigenvalues  Aj,  i  =  1,2, ,  and  their  corresponding  eigenfunctions 
(f)i,i  =  1,2, ... ,  can  be  determined  analytically  as  discussed  in  [36].  Different 
probability  distributions  can  be  chosen  for  1^.  The  effects  of  log  permeability 
with  uniform,  beta  and  Gaussian  distributions  on  the  mean  and  standard 
deviation  of  the  output  were  investigated  in  [37],  where  the  results  showed  that 
the  three  distributions  had  close  peak  values  of  standard  deviation.  Therefore, 
without  losing  the  main  feature  of  the  output  uncertainty,  here  17  are  assumed 
as  i.i.d.  uniform  random  variables  on  [—1,1]. 


In  this  problem,  the  fine-scale  permeability  is  defined  on  64  x  64  grid  and  the 
coarse  grid  is  taken  as  8  x  8.  For  comparison,  the  reference  solution  is  taken 
from  10®  MC  samples,  where  each  direct  problem  is  solved  using  the  fine-scale 
solver.  The  stochastic  problem  is  solved  using  HDMR,  where  the  solution  of 
each  deterministic  problem  at  the  collocation  points  is  from  the  multiscale 
solver.  In  this  way,  the  accuracy  of  both  multiscale  solver  and  HDMR  can  be 
verified.  In  our  previous  work  [43],  the  effects  of  the  correlation  length  and 
standard  deviation  have  been  studied  thoroughly.  Thus,  here  we  will  fix  the 
standard  deviation  to  cr^  =  1.0  and  investigate  the  effect  of  the  anisotropy  of 
the  random  field. 


5.3.1  Isotropic  random  field 

In  this  problem,  we  take  Li  =  L2  —  0.1.  Due  to  the  slow  decay  of  the  eigenval¬ 
ues,  Eq.  (70)  is  truncated  after  100  terms.  Therefore,  the  stochastic  dimension 
is  100.  The  problem  is  solved  with  HDMR  where  each  sub-problem  is  solved 
through  ASGC.  We  take  e  =  10“®,  =  5  x  10“®  and  6*2  =  10“^. 

In  Fig.  15,  we  compare  the  mean  and  standard  deviation  at  0.2  PVI.  It  is 
interesting  to  note  that  although  the  permeability  field  shows  heterogeneity 
for  different  realizations,  the  mean  saturation  is  the  same  as  the  solution  with 
homogeneous  mean  permeability  field.  This  behavior  is  called  “heterogeneity- 
induced  dispersion”  where  the  heterogeneity  smoothes  the  water  saturation 
profile  in  the  ensemble  sense.  Our  results  again  confirms  this  phenomenon, 
which  was  first  investigated  in  [31]  through  method  of  moment  equations. 
The  figure  also  indicates  that  higher  water  saturation  variations  are  concen¬ 
trated  near  displacement  fronts,  which  are  areas  of  steep  saturation  gradi¬ 
ents.  Therefore,  the  comparisons  between  the  MG  and  HDMR  results  are  only 
shown  around  the  displacement  fronts  on  the  bottom  two  plots  in  Fig.  15.  It  is 
seen  that  the  solutions  from  HDMR  compare  quite  well  with  the  Monte  Carlo 
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results.  The  convergence  of  HDMR  is  shown  in  Table  3,  where  the  normalized 
error  is  defined  the  same  as  before  with  MC  results  as  the  reference  solution. 
Ni  denotes  the  number  of  important  dimensions  and  N^.  denotes  the  total 
number  of  component  functions.  The  expansion  order  of  HDMR  for  all  three 
cases  is  2.  For  conventional  HDMR,  the  total  number  of  component  functions 
is  5051.  However,  by  using  adaptivity,  Nc  is  reduced  to  1047  which  clearly 
demonstrates  the  advantage  of  our  methods.  From  the  table,  it  is  seen  that 
the  results  are  indeed  quit  accurate  despite  the  fact  that  64-fold  upscaling 
is  used  to  solve  the  deterministic  problem  and  adaptive  methods  are  used  to 
solve  the  stochastic  problem. 


Fig.  15.  Mean  and  standard  deviation  of  saturation  at  0.2  PVI  for  isotropic  random 
field.  Top:  Mean  (a)  and  standard  deviation  (b)  form  HDMR.  Bottom:  Comparison 
of  mean  (c)  and  standard  deviation  (d)  between  MC  and  HDMR  near  the  satnration 
front. 

Table  3 

Convergence  of  HDMR  with  different  9i  at  0.2  PVI  for  isotropic  random  field. 


Oi 

Ni 

Nc 

#  Points 

Error  mean 

Error  std 

1  X  10“3 

2 

102 

1694 

7.47  X  10“^ 

4.38  X  10“2 

1  X  10"^ 

27 

452 

34379 

5.69  X  10-^ 

2.06  X  10-2 

5  X  10“® 

44 

1047 

77988 

5.10  X  10“^ 

6.66  X  10-3 
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Next,  we  demonstrate  the  interpolatory  properties  of  the  HDMR  method. 
As  mentioned  before,  one  of  the  advantages  of  HDMR  is  that  it  can  serve 
as  a  surrogate  model  for  the  original  problem.  Realization  of  the  saturation 
for  arbitrary  random  input  can  be  obtained  through  HDMR.  To  verify  this 
property,  we  randomly  generate  one  input  vector  and  reconstruct  the  result 
from  HDMR.  At  the  same  time,  we  run  a  deterministic  problem  with  the 
fine-scale  model  and  the  same  realization  of  the  random  input  vector.  The 
comparison  of  these  results  are  shown  in  Fig.  16.  In  addition,  in  Fig.  17,  we 
also  plot  the  probability  density  function  (PDF)  and  cumulative  distribution 
function  (CDF)  at  point  (0.2,0)  where  it  has  the  highest  standard  deviation 
as  indicated  from  Fig.  15(b).  These  results  indicate  that  the  corresponding 
HDMR  approximations  are  indeed  very  accurate.  Therefore,  we  can  obtain 
any  statistics  from  this  stochastic  reduced-order  model,  which  is  an  advantage 
of  the  current  method  over  the  MC  method. 


Fig.  16.  Prediction  of  the  saturation  profile  using  HDMR  and  the  solution  of  the 
deterministic  fine-scale  problem  with  the  same  input  for  isotropic  random  field. 
Left:  Saturation  at  0.2  PVI  from  direct  simulation  ,  Right:  Saturation  at  0.2  PVI 
reconstructed  from  HDMR. 


Fig.  17.  Isotropic  random  field:  (a)  PDF  of  the  saturation  at  point  (0.2,0)  and  0.2 
PVI,  (b)  CDF  of  the  saturation  at  point  (0.2,0)  and  0.2  PVI. 

Similar  results  at  0.4  PVI  are  also  given  in  Figs.  18,  19  and  20,  respectively.  It 
is  noted  that  the  standard  deviation  of  the  saturation  becomes  larger  at  later 
time  as  is  seen  from  the  wider  strip  of  the  non-zero  regions  in  the  contour 
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maps  at  0.4  PVI  in  Fig.  18.  With  the  increase  of  the  standard  deviation,  more 
collocation  points  are  needed  to  capture  the  overall  uncertainty.  Indeed,  there 
are  1229  component  functions  and  104662  collocation  points  in  this  case.  From 
Fig.  19,  it  is  seen  that  the  saturation  front  exhibits  a  much  more  signihcant 
variation  due  to  the  larger  standard  deviation.  Similarly,  in  Fig.  20,  we  plot  the 
PDF  and  CDF  at  point  (0.4,  0)  where  the  highest  standard  deviation  happens. 
It  is  noted  that  the  spread  of  the  PDF  at  0.4  PVI  is  wider  than  that  of  0.2 
PVI  which  again  indicates  the  larger  variation  of  the  saturation  at  this  time 
step.  Thus,  it  is  more  difficult  to  predict  the  uncertainty  with  the  simulation 
time  increases. 


Fig.  18.  Mean  and  standard  deviation  of  saturation  at  0.4  PVI  for  isotropic  random 
field.  Top:  Mean  (a)  and  standard  deviation  (b)  form  HDMR.  Bottom:  Comparison 
of  mean  (c)  and  standard  deviation  (d)  between  MC  and  HDMR  near  the  satnration 
front. 


5.3.2  Anisotropic  random  field 

In  this  problem,  we  take  Li  =  0.25,  L2  =  0.1.  Due  to  the  increase  of  the 
correlation  length  in  the  x  direction,  Eq.  (70)  is  truncated  after  50  terms. 
Therefore,  the  stochastic  dimension  is  taken  as  50. 

We  first  solve  this  problem  at  time  0.2  PVI  using  HDMR  with  ASGC.  We  take 
£  =  10“®,^!  =  5  X  10“^  and  62  =  10“^.  The  results  are  shown  in  Fig.  21.  It 


34 


0.9 

0.8 

0.7 

0.6 

0.5 


0.4 


I 


0.3 

0.2 

0.1 

0 


Fig.  19.  Prediction  of  the  saturation  profile  using  HDMR  and  the  solution  of  the 
deterministic  fine-scale  problem  with  the  same  input  for  isotropic  random  field. 
Left:  Saturation  at  0.4  PVI  from  direct  simulation  ,  Right:  Saturation  at  0.4  PVI 
reconstructed  from  HDMR. 


Fig.  20.  Isotropic  random  field:  (a)  PDF  of  the  saturation  at  point  (0.4,0)  and  0.4 
PVI,  (b)  CDF  of  the  saturation  at  point  (0.4,0)  and  0.4  PVI. 

is  interesting  to  note  that  the  shape  of  contours  is  nearly  the  same  as  that  of 
the  isotropic  random  field.  Only  the  values  of  standard  deviation  are  different. 
The  introduction  of  anisotropy  has  the  effect  of  increasing  the  output  uncer¬ 
tainty.  The  convergence  of  HDMR  shown  in  Table  4.  Again,  the  HDMR  re¬ 
sults  compare  very  well  with  the  reference  solution.  According  to  our  previous 
numerical  results  in  [43],  larger  uncertainty  requires  more  expansion  terms. 
Indeed,  more  expansion  terms  and  collocation  points  are  needed  compared 
with  that  of  isotropic  case.  In  addition,  the  highest  HDMR  expansion  order  is 
3.  There  are  3  third-order  component  functions,  which  indicates  the  existence 
of  higher-order  cooperative  effects  among  the  inputs.  The  reconstruction  of 
the  saturation  profile  is  shown  in  Fig.  22.  The  PDF  and  CDF  at  point  (0.2,  0) 
are  shown  in  Fig.  23. 


Finally,  we  show  that  HDMR  is  indeed  a  versatile  method  where  each  sub¬ 
problem  can  be  solved  by  any  stochastic  method.  Therefore,  we  solve  the 
problem  at  0.4  PVI  using  HDMR  where  each  sub- problem  is  solved  with  sparse 
grid  based  on  Gauss-Legendre  quadrature  rule  instead  of  ASGC.  A  level  3 
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Fig.  21.  Mean  and  standard  deviation  of  saturation  at  0.2  PVI  for  anisotropic  ran¬ 
dom  field.  Top:  Mean  (a)  and  standard  deviation  (b)  form  HDMR.  Bottom:  Com¬ 
parison  of  mean  (c)  and  standard  deviation  (d)  between  MC  and  HDMR  near  the 
saturation  front. 

Table  4 

Convergence  of  HDMR  with  different  9\  at  0.2  PVI  for  anisotropic  random  field. 


Oi 

Ni 

Nc 

#  Points 

Error  mean 

Error  std 

1  X  10“^ 

8 

79 

6199 

1.14  X  10“3 

4.69  X  10-2 

1  X  10"^ 

38 

754 

72243 

6.95  X  10"^ 

1.35  X  10-2 

5  X  10“® 

45 

1044 

96999 

6.51  X  10“^ 

1.01  X  10-2 

sparse  grid  is  chosen  for  each  sub-problem.  9i  is  chosen  as  1  x  10“^.  The 
results  are  shown  in  Fig.  24.  The  convergence  of  HDMR  is  given  in  Table  5. 
In  this  extreme  case,  all  the  50  dimensions  are  considered  as  important  and 
the  maximum  expansion  order  is  4.  This  again  is  consistent  with  our  previous 
results  in  [43].  Higher-order  terms  are  needed  to  capture  the  large  variability. 
Without  adaptivity,  there  are  251176  component  functions  for  a  4-th  order 
conventional  HDMR.  The  advantage  of  adaptive  HDMR  is  more  impressive 
in  this  case.  We  also  solve  this  problem  directly  with  a  50-dimensional  sparse 
grid  based  on  Gauss-Legendre  quadrature  rule.  The  results  from  levels  2  and 
3  sparse  grids  are  given  in  Fig.  25.  Since  the  mean  saturations  are  nearly  the 
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Fig.  22.  Prediction  of  the  saturation  profile  using  HDMR  and  the  solution  of  the 
deterministic  fine-scale  problem  with  the  same  input  for  anisotropic  random  field. 
Left:  Saturation  at  0.2  PVI  from  direct  simulation  ,  Right:  Saturation  at  0.2  PVI 
reconstructed  from  HDMR. 


Fig.  23.  Anisotropic  random  field:  (a)  PDF  of  the  saturation  at  point  (0.2,0)  and 
0.2  PVI,  (b)  CDF  of  the  saturation  at  point  (0.2, 0)  and  0.2  PVI. 

same,  we  only  show  the  comparison  between  standard  deviations.  For  level 
2  sparse  grid,  the  number  of  collocation  points  is  5301  with  the  mean  error 
8.31  X  10“^  and  std  error  4.38  x  10“^.  However,  when  increasing  the  sparse 
grid  to  level  3  with  a  total  number  of  192201  collocation  points,  the  mean 
error  increases  to  1.90  x  10“^  and  std  error  increases  to  7.09  x  10“^.  In  other 
words,  the  direct  sparse  grid  method  fails  to  converge.  It  is  computationally 
prohibiting  to  increase  the  sparse  grid  level  to  4  since  it  would  require  5402401 
collocation  points.  The  failure  of  convergence  is  due  to  the  steep  saturation 
gradient  near  the  displacement  front.  For  such  problems,  it  is  widely  known 
that  the  polynomial  based  quadrature  method  has  difficulty  in  convergence. 
From  the  results  shown,  it  seems  that  the  adaptive  HDMR  can  reduce  the 
irregularity  of  the  stochastic  space  through  decomposition  of  the  dimensions. 
However,  a  higher  order  expansion  may  be  needed  at  a  significant  increase  in 
the  computational  cost. 

Finally,  we  want  to  comment  on  the  computational  time  of  this  example. 
First,  in  Fig.  26,  the  convergence  of  standard  deviation  of  the  saturation  at 
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Fig.  24.  Mean  and  standard  deviation  of  saturation  at  0.4  PVI  for  anisotropic  ran¬ 
dom  field.  Top:  Mean  (a)  and  standard  deviation  (b)  form  HDMR.  Bottom:  Com¬ 
parison  of  mean  (c)  and  standard  deviation  (d)  between  MC  and  HDMR  near 
the  saturation  front.  Here  each  sub-problem  is  solved  using  sparse  grid  based  on 
Gauss-Legendre  quadrature  rule. 

Table  5 

Convergence  of  HDMR  with  different  d\  at  0.4  PVI  for  anisotropic  random  field. 


1 — 1 

Ni 

iVc 

Order 

#  Points 

Error  mean 

Error  std 

1  X  10“^ 

10 

96 

2 

4126 

1.32  X  10“3 

5.17  X  10“2 

1  X  10“^ 

38 

763 

3 

54925 

7.00  X  10“^ 

4.10  X  10“2 

5  X  10"® 

45 

1087 

3 

82407 

6.40  X  10"^ 

3.21  X  10-2 

1  X  10“® 

50 

2050 

4 

218136 

2.97  X  10“^ 

1.97  X  10-2 

one  point  with  the  number  of  MC  simulations  is  given.  The  points  are  chosen 
at  the  place  where  the  largest  standard  deviation  occurs  and  they  are  different 
for  different  cases.  From  the  figure,  it  is  seen  that  at  least  10^  MC  samples 
are  needed  in  order  to  achieve  statistical  convergence.  However,  there  are 
still  some  small  oscillations  after  it.  As  is  well  known,  the  MC  convergence 
rate  is  therefore,  to  ensure  a  good  comparison  with  HDMR,  we  use 

10®  samples  eventually.  It  took  about  19  hours  on  60  processors  while  the 
average  computational  time  for  HDMR  is  5  hours  on  the  same  number  of 
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Fig.  25.  Standard  deviation  of  saturation  at  0.4  PVI  for  anisotropic  random  field  us¬ 
ing  50-dimensional  sparse  grid  based  on  Gauss-Legendre  rule:  Comparison  of  stan¬ 
dard  deviation  between  MC  and  sparse  grid  level  2  (left)  and  3  (right)  near  the 
saturation  front. 

processors  in  such  a  high-dimensional  case.  It  is  also  noted  from  the  figure 
that  much  more  points  are  needed  to  achieve  statistical  convergence  in  the 
anisotropic  case  which  partially  explains  the  larger  variations  of  saturation  as 
was  seen  earlier.  Moreover,  an  interesting  observation  is  that  the  shapes  of  the 
convergence  plot  are  nearly  the  same  at  the  two  time  instants  for  the  same 
random  input.  This  phenomenon  suggests  that  although  the  convergence  rate 
of  MC  is  independent  of  the  number  of  stochastic  dimensions,  it  does  depend 
on  the  regularity  of  the  stochastic  input  space.  In  general,  more  MC  samples 
are  needed  for  a  stochastic  space  which  is  not  smooth  as  is  seen  from  the  case 
of  the  anisotropic  random  field. 


6  Conslusions 


In  the  first  part  of  this  paper,  a  new  multiscale  methodology  using  mixed  finite 
element  method  is  developed  for  the  solution  of  elliptic  equation  arising  from 
the  heterogeneous  porous  media  flow  problem.  This  multiscale  methodology  is 
based  on  the  framework  of  the  heterogeneous  multiscale  method  which  adds  a 
new  perspective  into  the  area  of  numerical  multiscale  methods.  A  novel  bound¬ 
ary  condition  for  the  local  cell  problem  is  proposed  which  gives  more  realistic 
flow  conditions  across  coarse-element  interface.  In  addition,  a  reconstruction 
method  for  the  fine-scale  velocity  is  also  proposed,  which  ensures  the  continuity 
of  the  mass  at  both  local  and  global  scales.  The  first  two  numerical  examples 
considered  verify  the  accuracy  of  the  new  method.  However,  as  a  first  step 
towards  this  new  method,  only  a  single-phase  flow  and  transport  problem  are 
considered.  Our  ongoing  research  includes  investigating  the  multi-phase  flow 
and  incorporating  the  multiscale  source  terms  and  well  modeling. 
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Fig.  26.  Standard  deviation  of  the  saturation  at  the  point,  where  the  largest  value 
occurs,  obtained  from  MC  simulations  versus  the  number  of  realizations. 


In  the  second  part  of  this  paper,  we  consider  the  uncertainty  quantification 
when  the  permeability  field  is  modeled  as  a  random  field.  The  newly  developed 
multiscale  method  is  used  as  a  direct  solver  within  the  framework  of  ASGC  and 
HDMR.  Our  numerical  results  in  Example  3  compare  well  with  the  MC  results 
with  fine-scale  solvers,  which  again  verifies  the  accuracy  of  both  multiscale  and 
HDMR  methods.  Our  study  confirms  the  interesting  phenomenon  that  the 
introduction  of  permeability  heterogeneity  leads  to  the  heterogeneity-induced 
dispersion.  The  obtained  results  also  indicate  that  the  HDMR  expansion  can 
serve  as  an  accurate  surrogate  model  for  the  underlying  stochastic  problem. 
Therefore,  our  ongoing  research  also  includes  using  this  stochastic  framework 
for  multiscale  permeability  estimation  as  an  extension  to  our  previous  work 
in  [61].  In  addition,  the  input  uncertainty  involved  in  the  current  work  is  only 
from  the  analytical  KL  expansion  with  known  covariance.  It  is  more  interesting 
to  consider  data-driven  stochastic  input  models  from  experimental  data  such 
as  in  [48]. 
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